# climdo-graphs.R
# generates graphs for the climdo paper
# Endre tvinnereim, Sept. 2016
library(stm)
# Clear all
rm(list=ls(all=TRUE))

load("stm_models.RData")
load("prepared_corpus.RData")
load("prepared_data.RData")

stmPrevFit <- stm_models[[5]]$runout[[4]]
labelTopics(stmPrevFit, topics=NULL, n=8)[2]

topicorder  <- c(3,5,6,1,2,4,7)
topiclabels <- c("Transportation", "Energy", "Attribution", "Reduction", 
                 "International", "Consumption", "Government")

covariates <- as.data.frame(stmPrevFit$setting$covariates$X)
colnames(covariates)

out$meta

str(stmPrevFit)

prep <- estimateEffect(1:7 ~ age+gender2+education2+education3+political_interest, 
                       stmPrevFit,
                       meta=stmPrevFit$setting$covariates$X)
# 2*2 box of charts
par(mfrow=c(2,2))
# Gender
table(covariates$gender2)
# genderplot <- 
stm:::plot.estimateEffect(prep,     main="Gender",
                    covariate = "gender2", 
                    topics = topicorder,
                    model=stmPrevFit, 
                    labeltype="custom", custom.labels=topiclabels, 
                    method="difference",
                    # verbose.labels=TRUE,
                    xlim=c(-.08,.08),
                    cov.value1=1, cov.value2=0, 
                    xlab="       Male        ...        Female", 
                    fg="grey")
# cov.value 1 is to the right; cov.value 2 to the left. See vignette p. 17. 

# Education
table(covariates$education2)
table(covariates$education3) # probably use this
table(out$meta$education)

# eduplot <- 
plot(prep,     main="Education",
                    covariate = "education3", 
                    topics = topicorder,
                    model=stmPrevFit, 
                    labeltype="custom", custom.labels=topiclabels, 
                    method="difference",
                    # verbose.labels=TRUE,
                    xlim=c(-.08,.08),
                    cov.value1=1, cov.value2=0,
                    fg="darkgrey",
                    xlab="       Not higher education       ...        Higher education")
# cov.value 1 is to the right; cov.value 2 to the left. See vignette p. 17. 

# Age
table(covariates$age)

ageplot <- stm:::plot.estimateEffect(prep,     main="Age",
                    covariate = "age", 
                    topics = topicorder[c(2,4,6)],
                    labeltype="custom", custom.labels=topiclabels[c(2,4,6)], 
                    method="continuous",
                    # verbose.labels=TRUE,
                    # xlim=c(-.1,.1),
                    model=stmPrevFit, 
                    fg="grey",
                    cov.value1=1, cov.value2=0, 
                    linecol=c("red", "black", "green"))
                    
                    # xlab="Age group")
# cov.value 1 is to the right; cov.value 2 to the left. See vignette p. 17. 

# Political interest
table(covariates$political_interest)

polintplot <- stm:::plot.estimateEffect(prep,     main="Political interest",
                    covariate = "political_interest", 
                    topics = topicorder[c(2,6,7)],
                    labeltype="custom", custom.labels=topiclabels[c(2,6,7)], 
                    method="continuous",
                    # verbose.labels=TRUE,
                    # xlim=c(-.1,.1),
                    model=stmPrevFit, 
                    cov.value1=1, cov.value2=0, 
                    linecol=c("red", "green", "blue"),
                    ylab="",
                    xlab="Political interest")
# cov.value 1 is to the right; cov.value 2 to the left. See vignette p. 17. 



plot(ageplot$x, ageplot$means[[2]], 
     main="Age",
     type="l", lty=1, ylim=c(.05,.25), xlab="", 
     ylab="Expected topic proportion", xaxt="n",
     fg="grey")
labels <- c("18-25","26-35", "36-45", "46-55", "56-65", "66-75", "75+  " )
text(1:7, par("usr")[1] - .685, srt = 90, adj = 2.5, labels = labels, xpd = TRUE)

lines(ageplot$x, ageplot$ci[[2]][1,], lty=2)
lines(ageplot$x, ageplot$ci[[2]][2,], lty=2)
text(3.5,.16, "Reduction")

lines(ageplot$x, ageplot$means[[1]], type="l", lty=1, col="red")
lines(ageplot$x, ageplot$ci[[1]][1,], lty=2, col="red")
lines(ageplot$x, ageplot$ci[[1]][2,], lty=2, col="red")
text(2,.245,"Energy", col="red")

lines(ageplot$x, ageplot$means[[3]], type="l", lty=1, col="darkgreen")
lines(ageplot$x, ageplot$ci[[3]][1,], lty=2, col="darkgreen")
lines(ageplot$x, ageplot$ci[[3]][2,], lty=2, col="darkgreen")
text(5,.06,"Consumption", col="darkgreen")

# Redone political interest plot
# 1: Energy (red)
plot(polintplot$x, polintplot$means[[1]], type="l", lty=1, ylim=c(.05,.25), 
     xlab="", main="Interested in politics", ylab="Expected topic proportion", xaxt="n", col="red",
     fg="grey")

labels <- c("Not at all","", "Somewhat", "", "Very" )
text(1:5, par("usr")[1] - .82,  adj = 0, labels = labels, xpd = TRUE)
lines(polintplot$x, polintplot$ci[[1]][1,], lty=2, col="red")
lines(polintplot$x, polintplot$ci[[1]][2,], lty=2, col="red")
text(2,.22,"Energy", col="red")
#2: Consumption (green)
lines(polintplot$x, polintplot$means[[2]], type="l", lty=1, col="darkgreen")
lines(polintplot$x, polintplot$ci[[2]][1,], lty=2, col="darkgreen")
lines(polintplot$x, polintplot$ci[[2]][2,], lty=2, col="darkgreen")
text(3,.06,"Consumption", col="darkgreen")
#3: Government (blue)
lines(polintplot$x, polintplot$means[[3]], type="l", lty=1, col="blue")
lines(polintplot$x, polintplot$ci[[3]][1,], lty=2, col="blue")
lines(polintplot$x, polintplot$ci[[3]][2,], lty=2, col="blue")
text(4,.16, "Government", col="blue")

